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We show that a coherent picture of the dc conductivity of monolayer and bilayer graphene at 
finite electronic densities emerges upon considering that strong short-range potentials are the main 
source of scattering in these two systems. The origin of the strong short-range potentials may lie 
in adsorbed hydrocarbons at the surface of graphene. The equivalence among results based on the 
partial-wave description of scattering, the Lippmann-Schwinger equation, and the T matrix approach 
is established. Scattering due to resonant impurities close to the neutrality point is investigated via 
a numerical computation of the Kubo formula using a kernel polynomial method. We find that 
relevant adsorbate species originate impurity bands in monolayer and bilayer graphene close to the 
Dirac point. In the midgap region, a plateau of minimum conductivity of about e 2 /h (per layer) is 
induced by the resonant disorder. In bilayer graphene, a large adsorbate concentration can develop 
an energy gap between midgap and high-energy states. As a consequence, the conductivity plateau 
is supressed near the edges and a "conductivity gap" takes place. Finally, a scattering formalism 
for electrons in biased bilayer graphene, taking into account the degeneracy of the spectrum, is 
developed and the dc conductivity of that system is studied. 



PACS numbers: 81.05.ue, 72.80.Vp, 78.67.Wj 
I. INTRODUCTION 

In his famous bookJUPeierls noted that in three dimen- 
sions the first Born approximation (FBA) cannot be used 
to deal with short-range potentials in general, even when 
the potential is not too strong. The reason lies in the 
fact that the FBA overestimates the value of the scat- 
tering cross section and modifies the energy dependence 
of the latter relative to the exact result. The fundamen- 
tal reason why this effect takes place has its roots in the 
modification of the wave function within the region where 
the potential is finite. There, even for moderate poten- 
tials, the wave function is strongly deformed relative to 
the plane wave used in the FBA. 

Since his main concern was nuclear physics, Peierls did 
not address the validity of the FBA in systems of reduced 
dimensions. Contrary to nuclear physics, some con- 
densed matter systems impose dimensional constraints 
on the electronic motion — a direct consequence of the 
lattice structure of the given solid. Electrons moving in 
graphene face the most dramatic dimensional constraint, 
being forced to move along a strictly two-dimensional 
plane formed by a honeycomb lattice of carbon atoms! 2 ^ 
In bilayer graphene, electrons are also confined to move 
in two dimensions. Since bilayer systems are a stacking 
of two graphene sheets, the electrons may, additionally, 
hop between the layers. 

Scattering cross sections in condensed matter physics 
are of ultimate importance for the interpretation of dc 
transport in solids, especially concerning the effect of 
localized impurities. These can be described by either 



short-range or long-range potentials. Following Peierls 
the correct interpretation of the conductivity of a metal 
at low temperatures may require a description of elec- 
tronic scattering by impurities beyond the FBA: this 
is particularly true if the impurities give rise to strong 
short-range potentials. 

In systems such as monolayer and bilayer graphene, 
where the electronic density can be tuned between and 
~ 10 14 cm _2 P computing the correct dependence of the 
cross section on the Fermi energy is a crucial ingredient 
for a meaningful interpretati on of the data. Since the 
early days of graphene physics it became clear that the 
conductivity of monolayer graphene shows slightly sub- 
linear dependence on electronic density. On the other 
hand, the conductivity of bilayer graphene shows, con- 
sistently, a robust linear dependence on the backgate 
potential. Both monolayer and bilayer graphene-based 
field-effect devices use sheets from flakes produced in ex- 
actly the same manner, i.e., via exfoliation of graphite. 
(More recently, graphene has been isolated via epitaxial 
growth onSiCP and chemical vapor depositions on metal 
surfaces P2I) It is now believed that the main sources 
of electronic scattering in exfoliated graphene are intro- 
duced during the device fabrication process. 

The sources of disorder in graphene can vary. They 
can be due to adsorbed chemical species, such as 
hydrogen atoms or hydrocarbon molecules, random 
strainpl rippling^L^I anc l scrolling,^ and electrostatic 
random potentials at the surface of the silicon oxide 
substrate caused by charged impurities EMUl (Chemi- 
cally synthesized graphene displays alternative scattering 
mechanisms P^J) 
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It is widely accepted that the strong carrier density 
fluctuations (electron-hole puddles) observed close to the 
neutrality or Dirac poi nt ar e due to localized subsur- 
face charged impurities! 15 * 19 ! Whether charged impuri- 
ties are also the limiting source of scattering in doped 
graphene (i.e. away from the neutrality point) remains 
unclear. In addition to charged scatterers, resonant scat- 
tering due to adsorbed hydrocarbons^ is currently as- 
cending in the list of candidates limiting the electr onic 
mobility in graphenepf^ll As we show in Sec. Ill A ad- 
sorbed hydrocarbons can effectively act as strong short- 
range scatterers. Strong, short-range, resonant scatter- 
ers can be mimicked by vacancies in a lattice modelPsHUl 
In magneto-optical transport studies of graphene, short- 
range scattering seems essential to explain the width of 
the cyclotron peak at high magnetic fieldPS 

Since the sources of scattering are likely introduced 
during the fabrication process, they must be the same 
for both monolayer and bilayer graphene. Therefore, a 
consistent theoretical description of the conductivity of 
graphene, at low temperatures and finite electronic den- 
sities, must be able to describe the experimental curves 
of both monolayer and bilayer graphene by invoking the 
same source of scattering. In this paper, we show that 
such a consistent theoretical description can be achieved 
by considering strong short-range potentials whose ori- 
gin may lie in adsorbed chemical species at the surface 
of the material. Instrumental to our description is the 
critical analysis developed by Peierls: Calculation of the 
exact scattering cross sections is essential for a correct 
interpretation of the experimental data. 

Before studying the dc conductivity for both mono- 
layer and bilayer graphene at finite electronic densities, 
a task we defer to Sec. |III[ we first survey the scatter- 
ing theory for electrons in these systems in Sec. [TT] This 
first step is essential for comprehension of the remaining 
sections. 

In Sec. |III| we show, using a simple and intuitive 
model, that the effect of adsorbed chemical species on 
graphene is equivalent to that of very strong on-site 
short-range potentials — the so-called resonant scatter- 
ers. Here, we use lattice-based numerical calculations of 
the density of states to show in some detail how this class 
of impurities affects the electronic structure of mono- 
layer and bilayer graphene. Using a continuous formula- 
tion, we also show that the semiclassical dc conductivity 
of both monolayer and bilayer graphene at finite densi- 
ties can be easily calculated using the intuitive approach 
to scattering given the partial-wave analysis. We apply 
the developed formalism to resonant scatterers, and show 
that this type of short-range disorder accounts well for 
experimental data. 

Further, we demonstrate the need for the computation 
of exact electronic scattering amplitudes when applying 
the Boltzmann approach to strong short-range potentials, 
an issue overlooked in the literature that we re-examine 
here. The validity of the semiclassical results at finite 
electronic densities and low impurity densities is estab- 



lished via a T-matrix calculation of the Kubo dc con- 
ductivity. Finally, by means of a numerical calculation 
based on the kernel polynomial method (KPM) , we illus- 
trate the breakdown of the semiclassical picture for elec- 
tronic densities close to the neutrality point. These sim- 
ulations explore the limit of finite impurity density, thus 
fully taking into account interference effects neglected in 
the Boltzmann approach. 

In Sec. |IV| we adapt the formalism of Sees. |IT]and |ITT| 
to describe scattering when a perpendicular electric field 
is applied to bilayer graphene. Conclusions are drawn in 
Sec. [V] Several technical aspects of our results are given 
in the Appendix. 

We note that transport in monolayer and bilayer 
graphene was addressed by some of us in an ealier 
publication. 55 However, it is important to remark that 
in the present work our goal is to provide a unified de- 
scription of transport in both systems based on the same 
scattering mechanism. Also, it is shown that the trans- 
port properties of the bilayer graphene can be understood 
in a much simpler, intuitive, and transparent way using 
the standard scattering formalism of partial waves. In 
this regard, our present work is complementary to the 
study developed in Ref. [55] That is, the present work 
closes the circle of showing that for both graphene and its 
bilayer, a coherent and unified description of dc transport 
in these systems can be described by one and the same 
formalism, be it the more formal and mathematically de- 
manding one of the transfer matrix or the intuitive and 
simple one of partial waves. 



II. PARTIAL- WAVE ANALYSIS IN GRAPHENE 

As discussed in Sec. [I] calculation of the dc conductiv- 
ity of a metal requires computing transport cross section 
as accurately as possible. A well-established approach 
is based on the computation of the phase shifts induced 
in the scattered electron wave function by the scattering 
potential. If the phase shifts are known exactly, so is the 
cross section. Below, we set the notation and introduce 
the central quantities needed in this work by giving a 
concise presentation of the phase-shift approach to scat- 
tering in the context of graphene and its bilayer J22H2U 
These results are later used in Sec. 



Ill Also, and to the 



best of our knowledge, the scattering theory for electrons 
in a biased graphene bilayer has not been developed so 
far in the literature, and therefore it is presented in Sec. 

CO 

Scattering theory states that the large-distance wave 
function of a particle in the presence of a scattering po- 
tential (with cylindrical symmetry) must have the form 
(in two dimensions) 



V>(r) 



/(*) 



(1) 



where ki = (fej,0) and kj = kt (cos 8, sin 9) are the mo- 
mentum of the incoming and scattered waves, respec- 
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tively; clearly, for elastic scattering, we have ki = kf = k. 
The scattering amplitude f(6) can be written in terms of 
the phase shifts 6 m associated with the partial-wave ex- 
pansion of the scattered wave function in the basis of 
angular momentum states. In Eq. (JT|), the first term 
represents the incoming particle, with the incoming mo- 
mentum oriented along the x axis, and the second one 
represents the cylindrical scattered wave function. 

As it stands, Eq. Q holds for the two-dimensional 
Schrodinger equation.^ However, for both monolayer 
and bilayer graphene, the large distance behavior of the 
wave function differs slightly, but significantly, from Eq. 
0. 



A. Electronic scattering in graphene 





Figure 1: (Color online) Lattice structure and Brillouin 
zone of monolayer graphene. Left: Hexagonal lattice of 
graphene, with the next nearest neighbor, Si, and the prim- 
itive, ai, vectors depicted. The area of the primitive cell is 
A c = 3\/3ao/2 ^ 5.1 A 2 , and a ~ 1.4 A. Right: Brillouin 
zone of graphene, with the Dirac points K and K' indicated. 
Close to these points, the dispersion of graphene is conical 
and the density of states is proportional to the absolute value 
of the energy. 

For graphene, the motion of the electrons in the 
7r orbitals is, at low energies, described by the two- 
dimensional massless Dirac Hamiltonian, reading^ 



H K = v F cr ■ p 



(2) 



where the Fermi velocity is defined as vf — 3tao/(2fi), 
t is the hopping integral between the p z orbitals of two 
adjacent carbon atoms, and a ~ 1.4 A is the carbon- 
carbon distance in graphene (see Fig. [T]). The vector <x 
is written in terms of Pauli's matrices as cr = {<j Xl a v ), 
and p is the momentum operator. The vector K denotes 
one of the two (inequivalent) edge points of the Brillouin 
zone, also called Dirac points or valleys. Because neu- 
tral graphene is half-filled (i.e., the ir orbitals contain 
one electron), these two points control the low-energy 
physics. Depending on the nature of disorder and the 
Fermi energy, coupling between momentum states from 
different valleys can take place. Intervalley scattering 
is known to induce weak localization corrections to the 
conductivity and, ultimately, fully localiz e states in the 
thermodynamic limit at zero temperaturefSSESEll 



In what follows, we assume that the two Dirac points, 
K and K', can be treated independently. This procedure 
is justified because intervalley scattering (known to occur 
for short-range scatterers) manifests itself primarily in 
the coherent regime, through backscattering interference. 
For low concentrations of scattering centers, finite sam- 
ple size, and finite temperatures (the typical experimen- 
tal situation), coupling between the Dirac points can be 
neglected when considering high enough electronic densi- 
ties. Hence, with the exception of the lattice calculations 
(Sees. HIE and III G I , we neglect intervalley scattering 



in the continuous model calculations and introduce the 
valley degeneracy index, g v = 2. In Cartesian coordi- 
nates, the eigenstate of the Hamiltonian in Eq. ^ has 
the explicit form 



*±(r) = 



1 



'2A 



1 

±e ie 



(3) 



with 8^ = arctan(fcj,/fc x ) and A denoting the total area of 
the system. The energy eigenvalues corresponding to the 
eigenfunction in Eq. ^ are E — ±Vphk. From the latter 
follows the density of states per spin and per unit cell, 
p{E) — 2\E\/ (it^/it 2 ), where the contribution from the 
two valleys has been taken into account. The probability 
density current reads^ 



(4) 



For the study of scattering, it is more convenient to recast 
the Hamiltonian in Eq. Q in cylindrical coordinates r 
and as 



Hk = —ivph 



L_ 
L+ 



(5) 



where the operators L± = e (d r ± ir 1 de) act as ris- 
ing/lowering operators, according to the following result 

L ± [C m (kr)e l9m } = TkC m±1 (kr)e^ m±1 ^ . (6) 

In Eq. ([6]), the function C m (kr) stands for J m {kr) and 
Y m (kr), the first-kind and second-kind Bessel functions, 
respectively, and for the Hankel functions of the first kind 
Hm^ and second kind Hm ■ For the modified Bessel func- 
tion K m (kr) we have 



L±[K m {kry 0m ] = -kK m±1 (kr)e 



■9(m±l) 



(7) 



In cylindric coordinates, the radial probability density 
current reads 



where a r is defined as 











(8) 



(9) 



The tangential component of the probability density cur- 
rent reads Jg — v F '$>\_<jg'fy ±, with ag = cr r diag(i, — i), 
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and where diag(i, — i) represents a diagonal matrix. Let 
us now derive, for massless Dirac electrons in two dimen- 
sions, the equivalent of the asymptotic wave function in 
Eq. ([lj. To that end, we note that a state having the 
form 



*m(r, 



1 



J m (kr) 
±ie ie J m+1 (kr) 



(10) 



is also an eigenstate of the Hamiltonian in Eq. ^j. We 
start by assuming that the asymptotic (large r) behavior 
of the wave function in the angular momentum channel 
m has the form (from here on, we consider only E > 0) 



1 



x e 



nAkr 



cos(kr - X m + S m ) 
ie sin(fcr — A m + S m ) 



(11) 



an ansatz inspired by the fact that the Dirac equation 
for graphene is a set of two coupled first-order differen- 
tial equations and in the asymptotic limit of the Bessel 
functions at large rP^ 



(12) 
(13) 



J m (x) = \ — cos(a; - A m ) , 

7TX 



-^m (^) 



sin (a; — A m ) 



with A m = m7r/2+7r/4. Using Eq. (Ill, we write the total 



wave function as an expansion in partial waves, reading 

oo 

*(r,0)= * m * m (r,e). (14) 

771 — — OO 

Exploiting of the relation 

oo 

= ]T i m e me J m (kr), 



Akr 



we obtain 



*(r) 



1 



1 



'2A 



1 



f(0) 



with the scattering amplitude reading 

f(o)-- 



J2 e m0 e l5 ™sm5 r , 



(15) 



(16) 



(17) 



It is a simple exercise to show that the first term in 
Eq. (16 1 corresponds to a flux J x = vp/A (and J y =0), 



whereas the second term corresponds to a radial flux 
J r = VF\f(0)\ 2 /(rA) (and Jg = 0). Thus, according to 
the usual definition of the differential cross section, a(9), 
it follows that 



<0) = |/(0)| s 



(18) 



Before we turn to scattering in bilayer graphene, it will 
be useful, for later use, to introduce other asymptotic 
forms of the Bessel functions J m (x), Y m {x), and K m (x), 
in addition to those already given in Eqs. (12 1 and (13 I. 
For large x, we have^ 



(19) 



K m (x) = 

For iCl, the asymptotic forms reacP^ 

J m (x) - (x/2) m T-\m+l), 

Y (x) = 27r~ 1 lna;, 



(20) 



(21) 



Y m (x) = -7r- 1 r(m)(x/2)- m , m = 1, 2, . . . , (22) 

and 

K (x) = -In a;, (23) 

K m (x) = 2- 1 r(m)(x/2)- m , m = l,2,..., (24) 

where T(x) is the gamma function. We now consider 
scattering in bilayer graphene. 

B. Electronic scattering in bilayer graphene 

Bilayer graphene has four atoms per unit cell, with the 
two honeycomb sheets arranged according to a Bcrnal 
stacking, as shown in Fig. [2] Two of the atoms belonging 
to each of the layers are on top of each other (atoms Ai 
and B2, in Fig. [2]), allowing for interlayer hopping. This 
pr ocess is represented by a hopping parameter, tj_ « 0.5 
e V! 36 l 37 l The other two carbon atoms, labeled A 2 and B 1 
in Fig. [2] are not coupled to the carbon atoms of the 
other layer, in accordance with the assumptions of the 
minimal model for electronic motion in bilayer graphene. 

The band structure of bilayer graphene has four bands, 
but the low-energy physics (\E\ <C t±) can be described 
by an effective model of only two bandsj2^H2H where the 
atoms linked by t± are projected out since they describe 
high-energy bands: the dimmer of atoms A\ and B2, 
linked by t±, form a two- level system with energy levels 
±t±. Additionally, the atoms in the two sheets can be 
made noncquivalent by applying an electric field perpen- 
dicular to the sheets, in this way inducing a gap in the 
spectrum (the electrostatic potential difference between 
the two layers is 2V)^ im 

The derivation of the effective Hamiltonian is straight- 
forward. We write the full Hamiltonian as 



H = 



V 








n 





-V 


TT-t 








7T 


-V 


-t± 


7ft 





-t± 


V 



rt 



H 



LH 



# LH #H 



, (25) 
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Figure 2: (Color online) Lattice structure of bilayer graphene. 
The atoms labeled Ai and B\ lie on the bottom graphene 
layer, whereas atoms A2 and B2 are in the top layer. Electrons 
can hop between layers via a perpendicular hopping parame- 
ter t± between the carbon atom A\ and carbon atom B2 (con- 
nected by solid lines). The Brillouin zone of bilayer graphene 
is the same as that of monolayer graphene (see Fig. [TJ . 



where the columns in the Hamiltonian are labeled by 
the four atoms in the unit cell. In ascending order, this 
labeling is B±, A 2 , B 2 , and A%. The operator tc is defined 
as 7r = vf (Px +ip y )- The eigenproblem H\ip) 
can be written as 



\<P) 

\x) 



E 



\<f) 

\x) 



E\il>) 
(26) 



It follows from Eq. ( 26 I that 



#i» + H Lli (E - H^hUv) = E\<p) 



(27) 



and considering that t± ^> (V, \E\), we have Hbl\<p) 
E\tp), witlP 



H 



BL 



Va z 



V 
if 



TTTT' 







— 7T 7T 



1 



(n) 2 



To keep things simple, in what follows we consider the 
case V = 0; later we discuss the case 7^0. In cylindric 
coordinates, the Hamiltonian, Eq. (28 1, is written as 



BL 



y\h 2 
ti 



Vl 




(29) 



It is important to stress two differences between the 
Hamiltonians in Eqs. ([2| and (29 1 regarding boundary 
conditions and the nature of the scattering states. To be 
concrete, let us assume that the electron is subjected to 
a potential well of the form V(r) = V$9(R — r). In the 
case of the Dirac Hamiltonian, the boundary conditions 
at r — R imply continuity of the two components of the 
spinors, whereas for the bilayer Hamiltonian we have to 
impose continuity of both the components of the spinors 
and their first derivative. The second aspect is related to 
the fact that elastic scattering conserves energy. Thus, 
since in bilayer graphene we have E = ±.v 2 F h 2 k 2 /t±, and 
keeping the energy constant, say E > 0, as in any scat- 
tering process, there are two admissible solutions: a real 
solution, k = \/t±E/(vph), and a purely imaginary one, 
k = iy/t±E I (vfK). Therefore, bilayer graphene supports 
evanescent modes at the interface r — R. This fact is es- 
sential to satisfy the boundary conditions obeyed by the 
wave function P^l 

As in the case of the Dirac Hamiltonian, we have to de- 
rive the form of the probability density current for elec- 
trons described by the Hamiltonian in Eq. (29). The 
usual procedure 34 gives that any component Ji of the 
current has the form 



where for I = x, y we have 



J x — @x &x 



Uydy , 



(31) 



(32) 



and 



J 1 1 



CTyd X 



'H "y*sx 0~ x dy . (33) 

For the radial component, t — r, we have 





Jr 



( 

,2i0/Q „-~-l; 



e M0 {d r - ir~ Y d e ) 

and for the tangential component, £ — 9, we have 

-ie- 2i9 {d r -ir- x d e ) 
ie 2i0 {d r + ir- l d e ) 



(34) 



J e = 



(35) 



Taking into account that the Hamiltonian in Eq. ( 29 1 



forms a set of two coupled second-order differential equa- 
tions, we assume that the asymptotic (large r) behavior 
of the wave function in the angular momentum channel 
m has the form 



and the eigenfunctions (regular at the origin) can be writ- 
ten as 



*m(r, 9) 



1 



>2A 



J m {kr) 
Te 2ie J m +2(kr) 



Am9 



(30) 



to which the eigenvalues E — ±v 2 F h 2 k 2 /t± correspond. 
From the latter result follows the density of states per 
spin and per unit cell, p(E) = tj_/(ny/3t 2 ), where we 
have included a factor of 2 coming from the valley 
degeneracvlTTI 



1 

nAkr 



cos(fcr - A m + 8 m ) 
e 2ie cos(fcr — A m + <5„, 



xe 



(36) 



Following the same procedure used to derive Eq. ([16| , we 
can show that the large-r behavior of the total electronic 
wave function in graphene bilayer in the presence of a 
potential has the form 



*(r) 



1 



'2A 



1 



'I A V e 



1 

JliQ 



ikr 



(37) 
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Using Eq. (31 1, we can easily conclude that the first term 



in Eq. (|37| corresponds to a flux J x = 2vphk/(At±) = 
v/A, where v is the velocity of the particle, and that 
the second term corresponds to a radial flux of the form 
J r = 2v%hk\f{6)\*/(rAt 1 _) = v\f(9)\ 2 /(Ar), with f{9) 
still given by Eq. (17 1. As before, it follows that the 
differential cross section is given by Eq. ( |T8| ). 

In Sec. |III| we apply the this formalism to the case 



of a potential well described by the potential V(r) — 
Vq6(R — r) in the strong interacting regime Vb 3> t. We 
will see that the results are insensitive to the particular 
form adopted for V(r) as long it corresponds to a strong 
short-range potential. 



III. THE DC CONDUCTIVITY OF GRAPHENE 
AND ITS BILAYER 

As discussed in Sec. |T] there is growing evidence that 
the limiting scattering mechanism of the electronic mobil- 
ity in graphene is due to strong short-range potentials, 
likely to have originated from adsorbed hydrocarbons. 
These adsorbed atoms and/or molecul es act as resonant 
scatterers, giving rise to midgap states .UnBDSSl 

This section is most important: it clarifies why the 
statement that short-range scatterers in graphene give a 
dc conductivity independent of the gate voltage is erro- 
neous. As noted in Sec. |T] this misleading idea has its 
roots in the FBA, which fails blatantly in this problem, 
as we demonstrate in what follows. 



A. Adsorbed atoms in graphene as strong 
short-range scattering centers 

The resonant scattering mechanism is easy to seize 
by considering a simple model. Let us write the tight- 
binding Hamiltonian of the n electrons in graphene as 
(spin index omitted) 



H=-tJ2\A,R n )(R n + S i ,B\+R.c. 



(38) 



where \A,R n ) represents the Wannier state at the unit 
cell _R„; an equivalent definition holds for \B,R n + Si), 
where 8i is one of three nearest-neighbor vectors in the 
honeycomb lattice, as depicted in Fig. [T] 

We now consider that an impurity is binding covalently 
to a carbon atom at site R n = 0. This situation adds to 



the Hamiltonian in Eq. ( 38 I a term of the form 



H IS = 04 d |ad)(A0| +h.c.) + e ad |ad)(ad| , (39) 

where V a( \ is the hybridization between the adatom (or a 
carbon atom of a hydrocarbon molecule) and a given car- 
bon atom of graphene, e a( j is the relative (to graphene 's 
carbon atoms) on-site energy of the electron in the 
adatom, and |ad) is the ket representing the state of the 



electron in the adatom. Taking the wave function to be 
of the form 

|V) = Y^l A (^)\ A > R n)+ B ( R n + S ^\ B > R n + d ^ 

n 

+ C ad | ad) , (40) 
the Schrodinger equation applied to the site R n — reads 
EA(0) - KdC ad = -t[B(6 1 ) + B(S 2 ) + B(6 3 )],(il) 



{E - e ad )C ad = Kd^(0). 
Solving for C a d, we obtain 



t[B(6!) + B{6 2 ) + B(8 3 )] = [E 



(42) 



ad 1 A(0) . 



E — e a d 



(43) 

The resonant effect is included in the last term in 



Eq. ( 43 1 , which represents an effective local potential of 



strength 



V^/{E 



£ad) 



(44) 



Quantum chemical calculation s can d etermine the value 
of the parameters e a d and Vad! 20 * 41 * ^ Typical values are 

U ad ~ It ~ 5 eV and e ad 0.2 eV, 1 ^ leading to 

V c s ~ 100 eV at half-filling (E = 0), a rather strong on- 
site potential. On the basis of this fact, it is natural to 
expect that adsorbates (i.e., resonant scatterers) and va- 
cancies lead to similar effects on the electronic structure 
and transport properties. In monolayer graphene, vacan- 
cies are known to significantly alter the density of states 
at energies close to e a d- In particular, vacancies induce a 
large spectral transference from the Van Hove singulari- 
ties to the neighborhood of the Dirac point. As a conse- 
quence, the densit y of s tates displays sharp peaks within 
the midgap region jSsEU This effect was first demonstrated 
in Ref. 26; recently, it has been sho wn th at indeed adsor- 
bates do originate similar behavior! 20 ! 72 ! 

Here we report similar results for bilayer graphene. 
To calculate the density of states, we employ the KPM 
(see Ref. for a review). For the sake of simplicity, 
we have considered equal concentrations of adsorbates in 
both bottom and top layers. (The actual applicability 
of this choice depends on the laboratory conditions and 
specific experimental setup.) In what follows, we discuss 
the situation where the adatoms bind only to carbons 
with coordination number z = 3 (i.e., those termed A 2 
and Bi in FigM. 

The effect of resonant impurities in the electronic 
structure of monolayer and bilayer graphene for differ- 
ent adsorbates concentrations, n a d, per carbon atom, is 
shown in Fig. [3] For illustration purposes, we present 
the results for a high defect concentration, n a d ~ 1%, so 
that the modification of the graphene electronic struc- 
ture is visible to the eye in a wide energy window; later 
we will see that the estimated values for defect concen- 
tration, for typical experimental conditions, are actually 



far below these values (Sees. IIIC and HID) 
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1.5r 



monolayer 

1 . Mayer (RS bind to A2 and B 1) 

2. biLayer (RS bind to any carbon) 





Figure 3: (Color online) Effect of adatoms (resonant scat- 
terers) on the density of states (DOS) of monolayer graphene 
(top) and bilayer graphene (bottom). Calculation of the DOS 
was carried out in honeycomb lattices with N = 1000 x 1000 
carbon sites for different concentrations of adsorbed atoms 
(periodic boundary conditions and 10 realizations of disorder 
were taken). The tight-binding parameters read V a d = 2i, 
£ad = -0.0625*, and t ± = 0.2*. The DOS discloses a dislo- 
cation of spectral weight toward the midgap region, a phe- 
nomenon first reported for vacancies in monolayer graphene 
in Ref. [26] 



In both graphene systems, the adatoms lead to well- 
defined peaks close to zero energy, the so-called midgap 
region. As mentioned above, such enhancement of the 
density of states is accompanied by a decrease in spec- 
tral weight near the Van Hove singulari ties, a situation 
reminiscent of vacancy- induced disorder ! 26 * 44 ! The effec- 



tive potential [Eq. ( 44 1 ] , despite being very strong, is 



bounded, explaining the slight electron-hole asymmetry 
near the Dirac point. The resonant peaks are centered 
at negative energies because e a( j < 0. Increasing the 
impurity concentration brings more spectral weight to- 
ward the midgap region. In bilayer graphene, though, 
a curious phenomenon takes place: when the impurity 
concentration is large enough, a gap opens separating 
midgap states, forming the impurity band, from high en- 



Figure 4: (Color online) Density of states (DOS) for bilayer 
systems with 5% resonant scatterers (RS) in the two scenarios 
described in the text, namely, (1) adsorbates binding only to 
carbons A2 and Bi, and (2) adsorbates forming bonds with 
carbons in any sublattice. The first situation opens a gap 
between the impurity band and high-energy states. The DOS 
of monolayer graphene is shown for comparison; tight-binding 
parameters are given in the caption to Fig. [3] 



ergy states (see Fig. [3] bottom) . Similar findings were re- 
ported in recent ab initio calculations considering asym- 
metric doping of graphene.^ 

Figure [4] shows how the electronic structure changes 
when the restriction on the allowable carbon-impurity 
bonds is relaxed. When adsorbates bind to carbons in 
any sublattice in the bilayer, the density of states is al- 
most indistinguishable from that of monolayer graphene 
(with the same impurity concentration) . The latter is ac- 
curate for a large energy window around the Dirac point 
(|e| < 0.5t); for higher electronic energies, the density of 
states becomes insensitive to the type of impurity-carbon 
bonds present in the bilayer samples. Roughly speaking, 
forming chemical bonds to every type of carbon decou- 
ples the layers, and hence dc-transport properties will be 
similar to those of a single layer of graphene (Sec. IIIG I. 

In light o f the present results and previous reports 
for vacancies 26 44 and resonant impurities in monolayer 
graphenej 20 * 72 ! we are led to conclude that the formation 
of an impurity band in the midgap region is universal 
in graphene systems with typical adsorbed species. In 
Sec. |III G| it will be shown that such an impurity band 
has a strong impact on the transport properties of un- 
doped graphene. 

Away from neutrality, the calculation of transport 
properties for the effective local poten tial mod el can be 
performed using the T matrix approach! 41 l 45 l 46 llts deriva- 
tion for resonant scatterers is elementary. It is well known 
that the T matrix for a local potential of intensity vq 
readfPHU 



T(E)=v [l-v G R (E)]- 



(45) 
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and 



G R {E) = ED~ 2 \n(E 2 /D 2 ) - nt\E\/D 2 



(46) 



with D ~ 3t. Then, using Eq. (44 1, the T-matrix due 



to an adatom must be of the form 



T(E) 



1 - V cS G R (E) E-e ad - V 2 d G R (E) 



(47) 



Since we are considering that V a d 2> (e a d, \E\), we can 
approximate the T matrix, Eq. (47 1, by 



T(E) 



1 



(48) 



G R (E) ' 

which is nothing but the T matrix for vacancies.^ 

The transport relaxation time r(kp) (at the Fermi sur- 
face) can be calculated using Fermi's golden rule, 



h/T{k F ) = Kn1\T(e F )\ 2 p{e F ) 



(49) 



where is the concentration of impurities per unit cell, 
and hp and ep are the Fermi momentum and energy, 
respectively. From the knowledge of t{Uf), the conduc- 
tivity of graphene follows from Boltzmann's transport 
equation (see the following section).^ 



B. The Boltzmann approach to dc conductivity 
using partial-wave expansion 

The above analysis made transparent that the effect 
of resonant scatterers is equivalent to that of a strong 
on-site potential (as long as the T-matrix formalism is 
applicable). We can then use the formalism of Sec. [TT] to 
compute the exact phase shifts in the presence of such 
a strong potential, from which r[kp) can be obtained. 
This type of calculations is equivalent, and alternative, 
to calculations based on the T-matrix approach in the 
lattice, with the appropriate choice of the effective size 
of the impurity. 

A relation between r(kp) and a(9) is provided bj0H 



l/r{k F ) = rii (v fcF • e r ) <t t , 



(50) 



where m is the concentration of impurities per unit area, 
Vfc F is the velocity of the electrons at the Fermi surface, 
e r is the radial versor in cylindric coordinates, and gj- is 
the total transport cross section,^ 



/•27T 

cr T = / d6 (1- cos 6)a(6) 
Jo 



(51) 



- J2 sm 2 (5 m -6 m+1 ) = -A(k). (52) 



m— — oo 



The conductivity of a given material follows from Boltz- 
mann's transport equation. The electric current has the 
general form 



(2tt)2 



dkr (fc)^M(v fc .E)v fc , (53) 



where n F is the Fermi distribution function, Sk is the 
dispersion of the electron, is the velocity of the particle 
with momentum k, E is the external electric field, and g s 
and g v are the spin and valley degeneracies, respectively. 
The electron velocity at the Fermi surface reads 



v F e r , 



whereas in the bilayer it has the form 



v fcF 



24 
t± 



hk F e r 



(54) 



(55) 



which depends on the position of the Fermi energy; the 
quantity M~ x = 2vp/t± plays the role of the electron's 
band mass. The dc-conductivity, 0^0 can be obtained 
from Ohm's law, j x — <j<i c E x . Combining Eqs. (52 1, (53 1, 
and ( 55 ) , the dc-conductivity for both monolayer 



and bilayer graphene has one and the same form, namely, 



Ode 



4c?- 



kp 



h 4n i A(k i 



(56) 



where the zero-temperature limit has been taken. The 
importance of Eq. (56) could not be more emphasized, 



since it shows that the final dependence of the conduc- 
tivity on k F , and therefore on the electronic density, is 
controlled by the behavior of A(kp), which depends only 
on the phase shifts S m ; these, in turn, depend on the na- 
ture of the scattering potential. Therefore, the exact cal- 
culation of the phase shifts emerges as the central theo- 
retical problem regarding the description of the variation 
of <7dc with the gate voltage for monolayer and bilayer 
graphene. 



C. Graphene 

For monolayer and bilayer, the electronic doping is con- 
trolled by a backgate voltage V g . The value of the Fermi 
momentum depends on the density of electrons, and, 
therefore, also on V g . If the dielectric between graphene 
(or its bilayer) and the backgate is made of silicon oxide 
and has a width of about 300 nm, then we have 



kp = iraVg , 



(57) 



■cm 2 ; numerically we have kp 



with a ~ 7.2 xlO 10 V" 
4.7 x 10- 3 x A" 1 . 

As we have discussed in Sec. |III A| an adsorbed atom 
or molecule (of specific types) can be described as an 
effective strong short-range potential. As a consequence, 
we model the effect of an adsorbed (resonant) chemical 
specie at the surface of graphene by a potential of the 
form 



V(r) = V 9(R ~ r) , 



(58) 



where R has to be of the order of ~lA and Vq ^> t. As a 
limiting behavior, we consider that Vq is made arbitrarily 
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large. In the Appendix we discuss the case where the 
potential is represented by a Dirac delta-function. The 
latter problem can be solved nonperturbatively, and an 
effective length scale i? e ff emerges in the problem due to 
an energy cutoff associated with the bandwidth. This 
effective length scale (R e g) is identified with the range R 
of the potential given above. Both problems lead to the 
same results for the conductivity of graphene (see later) . 

In the limit Vo —¥ oo, the potential defines an impene- 
trable barrier to the electronic probability flux. For elec- 
trons described either by the Schrodinger equation or by 
the Hamiltonian in Eq. ( 29 ) , the condition of zero flux 
for r < R is achieved by imposing that &(r = R) = 
[\]/(r) represents either a scalar or a spinor]. For electrons 
described by the massless Dirac equation, the latter im- 
plies that the wave function has to vanish everywhere 
and, therefore, cannot be used. In contrast, from Eq. ^ 
it is clear that the radial flux at r = R can be made 
if one of the components of the spinor is at r = 
In conclusion, the correct boundary condition enforcing 
zero flux at r = R for electrons in monolayer graphene is 
given by 



9i(r = R)=0, 



(59) 



where with i — 1,2, is one of the components of 
the spinor. Given the presence of two Dirac cones in 
graphene, it is immaterial which component we choose 
to obey the condition of Eq. 59 1 , as long as we consider 



the contributions to the two Dirac cones in the Brillouin 
zone of the honeycomb lattice. 



To satisfy the boundary condition in Eq. ( 59 I , we write 



the wave function describing the electrons being scattered 
by the barrier as 



e w J m+1 {kr) 



At 



Y m (kr) 
e ie Y m+1 (kr) 



Thus, the boundary condition in Eq. (59 1 implies that 

J m (kR) 



A? 
Af 



Y m (kR) ■ 



(60) 
at 

(61) 



Since for large r, the wave function in Eq. ( |60| must have 
the general form shown in Eq. (11), it follows that the 



ratio A™ I A™ has to be interpreted as 



■"-2 

A^ 



tan 8„ 



(62) 



which defines the phase shift S m . [A comment about the 
latter result is in order: In graphene, radially symmet- 
ric potentials originate phase shifts obeying S m = 8— m —i. 
This can be seen by noting that replacing m by — m— 1 in 



Eq. ( 10 1 produces another eigenstate of the Dirac Hamil- 



tonian. Equations (|61| and (62) show that impenetrable 
barriers force a different symmetry: S m — 5- rn .\ For 
backgate voltage values in the range V g < 100 V, and 
considering R ~ lA , we have Rk < 1 (known as the 
low-energy scattering regime). In this regime, the scat- 
tering is dominated by the s wave phase shift; that is, 
the dominant contribution to A(A:) comes from 



J (kR) 7T _ 

tan dn = — ) r- m — In 

° Y (kR) 2 



\kR), 



(63) 



1 i 1 i 1 i 



1/p. 




-45 -30 -15 15 30 45 
gate voltage, V 



40 -20 20 40 
gate voltage, V 



Figure 5: (Color online) Experimental data on graphene's 
conductivity. Left: raw data on a measurement of the resis- 
tivity, Pmeasured, of an exfoliated graphene sheet. Right: Fit 
of the conductivity, <7 S ub = 1/pmeasured, using Eq. |64|. The 



value of R was taken to be of the order of an and the fit pro 
vided an areal density of impurities of m ~ 2.5 x 10 11 cm"' 



where Eqs. (20 1 and (21 1 have been used. It follows from 
Eqs. (|52| and (63l that the conductivit y of graph ene ob- 



tained from Eq. 



has the final form 25 31 42 51 



4e 2 



Ode 



kp 



h 2tt 2 7 



\n {k F R). 



(64) 



Given that the value of R is constrained to be of the or- 
der of 1 A , n, is the only fitting parameter. Equation 
( 64 1 was used to fit the conductivity data^ of an exfoli- 
ated graphene sheet, as shown in Fig. [5] Because we took 
the limit Vo — > oo, the computed conductivity does not 
break electron-hole symmetry. The electron-hole asym- 
metry shown by the experimental data in Fig. [5] can be 
attributed to the presence of charge scatterers and / or to 
the role of the contacts.^! If we increase the value of R 
somewhat, the concentration of impurities needed to fit 
the data decreases. In Fig. [5] we have chosen to fit the 
conductivity for a positive gate voltage; it is manifest 
that Eq. (64) fits the data accurately [dashed (black) 
curve]. If we had decided to fit the data for negative 



(or a concentration n ad « 5 x 10 '' per carbon atom). In both values of V g , the obtained concentration of impurities 

panels, the horizontal dashed line (red) stands for twice the n ^ wou ld have been slightly different. The concentra- 

quantum of conductance, that is, 2e 2 /h. (Data from S. V. tion of sca tterers is rather small (see caption to Fig. IB] 

Morozov et al-H courtesy of A. K. Geim.) and agreeg ^ thg concentration of atomic scale defec V s 
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estimated via Raman measurements.^ This testifies to 
the strong effect of a few resonant scatterers dilluted in 
the surface of graphene (similar to atomic vacancies), as 
discussed in Sec. IIII Al 



The result given by Eq. (64) for the conductivity of 



monolayer graphene can also be obtained from a model 
where vacancies act as scattering centers.^ In view of 
the arguments given in Sec. |III A| this result comes as 
no surprise, since the effective local potential created 
by adsorbed hydrocarbons is much larger than the hop- 
ping integral t. Numerical simulations of the dc con- 
ductivity based on Kubo's formula in the presence of lo- 
cal potentials found a sublinear behavior for a graphene 
monolayer in qualitative agreement with Eq. (64). 



Let us now extend the previous analysis to the case of 
a graphene bilayer. 



D. Graphene bilayer 

Assuming that the dominant source of scattering in 
graphene is due to strong short-range potentials, then 
the same must be true for bilayer graphene. As a conse- 
quence, a consistent description of electronic scattering 
in both monolayer and bilayer graphene must use the 
same scattering potential to explain the measured con- 
ductivity in both systems. In the spirit of this work, this 
means that the scattering potential in Eq. ( 58 1 must also 



be used to compute the conductivity of graphene bilayer. 




V b (volt) 



Figure 6: (Color online) Dependence of the phase shifts 5o(-2) 
(solid line) and <Si(-3) (dashed line) on V g , for bilayer graphene 
with R — ao- The differences between the exact expressions in 
Eq. ( |66[ ) and the asymptotic values in Eqs. ( |69[ ) and ( |70[ ) are 
not visible to the eye. Other phase shifts are approximately 
within the same range of V g . 



As in the case of Eq. ( 60 ) , we seek a wave function in 



kinds, which in the present case assumes the form 



-4', 



Al 



Al 



-e 210 J m+2 {kr) 

Y m (kr) 
~e 2l0 Y m+2 {kr) 

K m (kr) 
~e 2l6 K m+2 {kr) 



(65) 



The introduction of the modified Bessel function K m (kr) 



in Eq. (81 1 is necessary to satisfy the boundary condition 
v]/(r = R)= 0. We recall that Hamiltonian in Eq. ((29 1 



supports evanescent waves at the boundary r = R, as 
discussed in Sec. II B Furthermore, for large r, K m (kr) 
decays exponentially, as we can see from Eq. ( 19 ). There- 



fore, at large distances, the behavior of the wave func- 
tion in Eq. (81 1 depends only on the form of J m (kr) 
and Y m (kr), as given by Eqs. (12 1 and (13). As a con- 
sequence, the phase shift 5 m is determined by the ratio 



A™ jA\ 



that is, we must have 
A 1 ? 



tan(L 



(66) 



as in the case of electrons in monolayer graphene [see 
Eq. (62 1], Imposing the boundary condition ty(r = R) = 



on the wave function (81 ), we obtain 

= ATJm(kR) + A™Y m {kR) + A™K m {kR) , (67) 
= ATJ m+2 (kR) + A™Y m+2 {kR) + A™K m+2 (kR) , 

from which follows 



A? = J m (kR)K m+2 (kR) - J m+2 (kR)K m (kR) 
A™ K m (kR)Y m+2 (kR) - K m+2 (kR)Y m (kR) 



(68) 



Combining Eqs. fl66| and (68 I, the equation for the phase 



shift S m follows at once. Contrary to the case of mono- 
layer graphene, the cross section is no longer dominated 
by 5q alone. The asymptotic expansions for 6q and 5i are 
(k F R < 1) 



tan<5 = - Mk F R/2)+j E - 1/2]- 



and 



tan S 1 = - [\n(k F R/2) + lE - l/4]-\ 



(69) 



(70) 



where 7^ = 0.577 ... is Euler's constant. In addition, we 
have two more nonzero phase shifts: 



5-2 = <5 , and <5_ 3 = Si. 



(71) 



These expressions are exact and reflect a symmetry of the 
eigenstates of Eq. ( 29 1 when radially symmetric scalar 



potentials are considered, namely, S m = 5- m - 2 - 

The dependence of 5q and 5\ on V g is given in Fig. [6| 
From Eqs. ([69]) (|7ll , it follows that A(k F ) ~ 4. The dc 



the form of a superposition of Bessel functions of different 



conductivity of bilayer graphene is, therefore, given by 

Ap 2 k 2 

aAc = ^-J^, (72) 
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E. T matrix approach for bilayer graphene 




-75 -50 -25 



25 50 75 



50 -25 
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Figure 7: (Color online) Fit of the conductivity data of bilayer 
graphene [solid (red) curve] using Eqs. \72\ and ( |75[ ). The fit 
has only a single parameter, the concentration of impurities. 
The obtained value is m ~ 4 x 10 10 cm" 



concentration 



Tlad : 



1 x 10 per carbon atom), for Eq. (72 1, and m ~ 1 x 



10 11 cm 2 (concentration n a( j ~ 0.25 x 10 5 per carbon atom) 
for the T matrix approach, using a model of pure vacancies, 



Eq. (751. Left: Data taken at a temperature of 20 K. Right: 



Conductivity of the same sample at the higher temperature 
of 100 K. The position of the Dirac point, Vd, was shifted to 
in this figure. (Data from S. V. Morison et al.r^ 
of A. K. Geim.) 



courtesy 



by Eq. ( 72 1. We have already shown that the effect of 



resonant scatterers can be captured by a model of pure 
vacancies, using both the T matrix and the partial-wave 
approaches. We now revisit the T matrix calculation in 
bilayer graphene^^ and show that, as in the case of the 
monolayer, a model of pure vacancies in the bilayer also 
captures the physics of resonant scatterers. 



In Refs. EH and |55J the calculation of the dc conduc- 
tivity took into account the full band structure of the 
graphene bilayer. That calculation could distinguish the 
four carbon atoms in the unit cell. In this section, we as- 
sume that vacancies are located at the two carbons that 
are not coupled by t±. 

In the notation in Refs. [Ml [55J the zero-temperature 
dc conductivity obtained from Kubo's formula is given 

by 



8c 2 
nh 



d(k 2 ) {lm[g» A (E F , k)]Im[g» B (E F + 5, k)} 



ImLO^F, k)]Sm\sfi%(E F + 5, k)}} 



(73) 



in the limit 5 — > 0; see Ref. 55 for the definitions 
of the Green's functions g(E,k). The fc 2 -integral can 
be performed exactly, as explained in Appendix C in 
Ref. [55J The resulting complicated formula can be ap- 
proximated by going through the following steps: (i) ne- 
glect the real part of the self-energies, (ii) expand the 
result in powers of the imaginary part of the self-energies 
r a (e) = — Im[E (e)], and (iii) assume that the energies 
involved fulfill \^\,tj_ ± > T A (e), T B (e). The leading 
term in this expansion yields the approximate formula 



2e 2 



Ode 



Ep(Ep + £j_) 



h t ± T B (E F ) + E F [T A (E F ) + T B (E F )Y 



(74) 



Curiously, the symmetry of the scattering amplitudes 
combine to make A.(k F ) independent of k F (with an accu- 
racy better than 1% in the relevant range of k F and R), 
making the conductivity proportional to the gate voltage. 
This result, together with the constant density of states 
(valid when \E\ <C t±), is at the heart of the exact linear 
dependence of the conductivity on the gate-voltage. We 
have used Eq. (72 1 to fit the conductivity data of an ex- 
foliated bilayer graphene sample, as shown in Fig. [7] The 
fit provides a concentration of impurities of the order of 
m fa 4 x 10 10 cm~ 2 (i.e., a concentration of adatoms per 
carbon atom of about n a d ~ 1 x 1CP 5 ). Since in bilayer 
graphene only two of the four surfaces are exposed to the 
environment, the rii value found above, being slightly 
smaller than that found for monolayer graphene, can be 
interpreted as a manifestation of this fact. 

Within the T matrix approach, the dc conduct ivity 
of bilayer graphene has been computed in the past!^^ 
The impurity concentrations used in those works were 
far too large to reveal the linear behavior in V g given cm" 



This expression is a good approximation for low impu- 
rity concentrations and away from the neutrality point, 
where the condition in step iii breaks down. This re- 
sult may be further simplified using the relation be- 
tween the Fermi energy and the density (assuming 
n, E F > 0) coming from the dispersion relation E F = 
\J (ij_/2) 2 + n(hv F ) 2 n — t±/2, resulting in 



2e 2 



n(hv F ) 



h t x T B (E F ) + E F [T A (E F ) + T B (E F )} ' 



(75) 



where n is the electronic density. To the extent that the 
denominator is independent of E F , the conductivity is 
linear in the density of carriers, n, in agreement with 
the description based on the phase shifts. For low im- 
purity densities, as is the case in exfoliated samples, the 
difference between the conductivity obtained from the co- 
herent potential approximation and the T matrix is very 
small except in a tiny region near the neutrality point. 
Using Eqs. (72) and (75 1, the data in Fig. ^ can be rea- 
sonably fit considering a density of vacancies of n., ~ 10 11 



F. Exact amplitudes versus first Born 
approximation 

The use of the FBA within the semiclassical Boltzmann 
approach is a common practice in condensed matter. In 
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dc-conductivity 


monolayer 


bilayer 


S 

FBA 


Const. 


~ k% 


nonperturbative; largeVo 


- [k F ln(k F R)] 2 


~ k% 


Hard-disk radius 7? " 


- [k F ln{k F R)] 2 





Table I: The conductivity due to a delta 8 potential: the 
FBA and the nonperturbative result in the relevant regime 
Vb 2> all scales. For comparison, the hard-disk result is listed. 
Although for the bilayer both the FBA and the exact calcu- 
lation give a conductivity proportional to k F , we should note 
that in the former case the conductivity is proportional to the 
strength of the potential, and therefore the FBA cannot be 
trusted in the regime of strong potentials, and the agreement 
of the two approaches is fortuitous. 



the present context, the FBA has been employed to inves- 
tigate the inter play between short-range and long-range 
scattering ! 57 ! 58 ! Jt s U se, however, requires the weak scat- 
tering condition to be verified. We have seen in Sec. |III A| 
that adsorbed atoms in graphene give rise to strong local 
potentials Vp > t , rendering inappropriate the use of 
the FBA for a description of scattering due to realistic 
short-range potentials. 



The form of the graphene conductivity [see Eq. ( 64 1] 
is not peculiar when hard-wall boundary conditions are 
present; potentials characterized by delta functions in 
real space yield equivalent results if exact scattering am- 
plitudes are considered instead of the FBA (see Ap- 
pendix). Moreover, beyond Boltzmann's kinetic the- 
ory, tight-binding calculations for graphene sheets with 
~ 0.02 /jm 2 show quantitative agreement with Eq. (64 1, 



while at the same time displaying qualitative disagree- 
ment with the FBAP 

To demonstrate that S potentials also mimic the ef- 
fect of strong range potentials, we calculate the exact 
scattering cross sections using the Lippmann-Schwingcr 
equation, an approach well suited to 5 potentials. (To 
the best of our knowledge, the case of bilayer graphene 
has not been considered before.) The calculations are 
shown in the Appendix and important limiting cases are 
summarized in Table HI 

For monolayer graphene the conductivity due to a S 
potential with strength Vq reads 



Ode 



4e^J2_ 
h rii 



{k F /2Ti)\n{k F R) 



(76) 



where R is a length scale introduced to regularize the 
Green's function. The FBA is recovered from Eq. (76) 



by considering Vq smaller than relevant scales, yielding 
a conductivity that does not depend on the carrier den- 
sity/gate voltage. In contrast, the strong scattering limit 
Vo 3> \E\ gives the same dependence found for the hard 
disk model [Eq. ( 64 1] upon the identification of R with 
the potential range. 

The situation is quite different in bilayer graphene be- 
ing described by a low-energy theory of massive electrons: 



both weak and strong scattering regimes yield a con- 
ductivity proportional to kp in the entire carrier density 
range; the exact result reads 



Ode 



4e 2 



1 



h 16n, 



8v 2 F h 2 
V t x 



kp. 



(77) 



Although for the bilayer both the FBA and the exact 
calculation results in the conductivity being proportional 
to kp, we should note that in the former case the con- 
ductivity is proportional to the strength of the poten- 
tial, and therefore the FBA cannot be trusted in the 
regime of strong potentials, and the agreement on the k in- 
dependence of the two approaches is fortuitous. (We re- 
mark that the limitations of the FBA for a description of 
electronic transport are not exclusive to short-range scat- 
terers and can also be found in Coulomb scatterersP^) 

The results of the present and previous sections con- 
firm the intuitive idea that delta potentials and hardwall 
(hard-disk) boundary conditions originate the same de- 
pendence of Ode on the Fermi momentum. Remarkably, 
letting Vq — > oo in Eqs. (76 1 and (77), give precisely 



Eqs. (64) and (72 1, respectively, and hence the two mod- 
els are equivalent with regard to strong short-range po- 
tentials. 



G. Quantum corrections near the neutrality point 

The Boltzmann approach beyond the FBA provides 
a good description of the effect of strong short-range 
scatterers on the transport properties of graphene at fi- 
nite carrier densities (and for not too large concentra- 
tions of resonant impurities) However, near the neu- 
trality point quantum interference effects become impor- 
tant and a fully quantum calculation is needed to as- 
sess dc-transport. (For recent reviews on the impor- 
tance of quantum effects in the transport properties of 
graphene see Refs. |5J [BUJ) In what follows, we present 
large-lattice, tight-binding numerical calculations in the 
low-density regime and finite (high) impurity concentra- 
tion limit n a d ~ 1%, where quantum corrections due to 
multi-scattering events cannot be ignored. 

Monolayer graphene — We start by extending the 
monolayer tight-binding Hamiltonian [Eqs. p8| and p9| ] 
to include a finite number AT a( j of adsorbed atoms of the 
same species, binding to carbons placed at (random) po- 
sitions {sj (i = 1, ...,N ad ), 

Hth = -t'£ i \R n ,A)(R n + 6i,B\+-R.c. 

n,Si 



^[Kdls^adXs^ai+H, 



2 = 1 



+e ad |s i ,ad)(s l ,ad|] 



(78) 



where Ci — A(B) for adatoms binding to carbon atoms 
in the A(B) sublattice. The Kubo formula for the zero- 
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Figure 8: (Color online) Conductivity as function of the nor- 
malized carrier density n c /n a d for a monolayer honeycomb 
lattice with TV = 1000 x 1000 for different concentrations of 
adsorbed atoms (periodic boundary conditions and ten real- 
izations of disorder were taken). The tight-binding parame- 
ters read V a d = 2t and e a d = —0.06254 . 



temperature dc-conductivity tensor reads^ 



<?ab 



(E) 



2nhe 2 



A 



-Tr 



v a 6(E-H th )v b 5(E~H< 



ti.> J 



(79) 

where v a ^ is the a(6)th component of the velocity op- 
erator (denned through the Heisenberg equation of mo- 
tion for the position coordinate) and A stands for area of 
graphene. 

We evaluate the longitudinal component of the conduc- 
tivity a xx employing a KPM: details of the calculation are 
given elsewhere. 62 The KPM amounts to approximate 
functions defined in bounded intervals by a truncated 
sum over polynomials with optimized weig htsPTo illus- 
trate the change in the transport properties near the neu- 
trality point, we simulate mesoscopic-size square sheets 
of graphene with N = 10 6 carbon sites. An adequate 
polynomial expansion of Eq. ( 79 ) allows us to perform 



the simulations with modest computational resources. 

We found that the expansion of Eq. ( [79] ) in Chebyshev 
polynomials of the first kind converges for concentrations 
of resonant impurities, n ad = N a ^/N, above a critical 
value n* d of about 1% (for N = 10 6 ). We interpret this 
result as an indication that for n a( j < n* d , electronic 
carriers are in the ballistic regime. (Recall that only in 
diffusive or localized regimes can a thermodynamic con- 
ductivity be defined.) The values n a( j > n* d correspond 
to concentrations of short-range scatterers several orders 
of magnitude larger than what is found in typical labo- 
ratory environments (about 10 _3 %; see previous sections 
and Ref. I2.'j|| but can, in principle, be reached via hydro- 



genation of graphene on SiC>2 The critical value n* d 
likely indicates the onset of diffusive behavior, I < L, 
where I is the mean free path and L denotes the lattice 
linear size. Thus, in principle be lowered by increasing 



L. 

Figure [8] shows results for conductivity as function of 
the carrier density; the latter was obtained by integration 
of the density of states p(E) (shown in Fig. [3]) , according 
to 



n c (E F ) = (g./D) 



p(E)dE, 



(80) 



where D = N + N a a is the total dimension of the prob- 
lem. The most peculiar feature in Fig.[8]is the plateau of 
finite conductivity, due to the formation of a low-energy 
impurity band (Fig.[3j top) , a particular case of disorder- 
enhanced conductivity ! 65 ! 66 1 7 ^ 

The dc conductivity at the neutrality point differs sig- 
nificantly from calculations based on Boltzmann kinetic 
theory. (1) The conductivity saturates at a low carrier 
density to a finite value er m i n > around e 2 /h (the pre- 
cise value depends on n ad and sample size) , in accordance 
with theoretical predictions.^ The width of the satura- 
tion is roughly proportional to the density of adatoms in 
the probed range of impurity concentration n a( j < 5% (a 
similar behavior was first reported using a self-consistent 
approximation to the Green's function of the electrons in 
the presence of a strong disordered potential and re- 
cently reported in Ref. 20). (2) The conductivity (for 
a fixed carrier density or energy) is not proportional to 
1/^ad- [In fact, a careful inspection of the KPM conduc- 
tivity data discloses that the latter observation extends 
to higher carrier densities: (resonant) adsorbate-limited 
transport in small samples of graphene displays a rich 
behavior until full diffusive transport is reached.] Both 
fact 1 and fact 2 clearly indicate that we are operating 
outside the applicability of the Boltzmann approach. 

Our results, in general, agree well with those reported 
in Ref. 20 for larger lattices (where N of the order of 10 s 
was used). Notwithstanding, we point out some differ- 
ences concerning the plateau of conductivity minimum: 
we observe neither peaks within the conductivity plateau 
(including for rt ad = 5%) nor a plateau's width of 2 x n ad , 
as claimed in that work. This could be due to the differ- 
ent methods and system sizes used (although in simula- 
tions with a larger lattice, we found no evidence of both 
effects) . 

A comment about intervalley scattering in our simu- 
lations is in order; Anderson localization induced by in- 
tervalley scattering will become experimentally relevant 
and prevent conductivity saturation only for either very 
strong disorder (i.e., high defect densities) or exceedingly 
large samples at very low temperatures. In contrast, our 
results, and those in Ref. [2TJJ show no evidence for local- 
ization even for relatively high amounts of resonant dis- 
order. This suggests that the localization length due to 
resonant scatterers is far larger than that obtained for an 
on-site Anderson model, hence allowing for conductivity- 
induced disorder, ao > 0, in typical-size graphene sam- 
ples. 

Bilayer graphene - The tight-binding Hamiltonian for 
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Figure 9: (Color online) Top: Conductivity as a function of 
the normalized carrier density n c /n a d for a bilayer honeycomb 
lattice with N — 2x 1000 x 1000 for different concentrations of 
adsorbed atoms (periodic boundary conditions and 10 realiza- 
tions of disorder were taken). The tight-binding parameters 
read Kd = 2t, e ad = -0.0625t, and t± = 0.2t. Bottom: 
Comparison of the conductivity (per layer) of monolayer and 
bilayer graphene (with n a d = 5% in both cases). Two bilayer 
curves are shown corresponding to different arrangements of 
resonant scatterers (RS) as discussed in Sec. ( |III A[ ): (1) adsor- 
bates binding only to carbons A2 and B\ , and (2) adsorbates 
forming bonds with carbons in any sublattice. The former 
situation leads to a supression of the plateau near the edges. 



bilayer graphene with resonant impurities reads 



-^tb LG) - #tb 1,2) +t±y^ j (\R n ,Ai)(R n ,B 2 \ 



,(L=1,2) 
'tb 



H.c.) 



n, Si 



J2 [Kd|sf,ad)(sf,C L |+H.c. 



i=l 
(£ = 1,2) 



-e ad |sf ,ad)(sf ,ad|] , 



(81) 



where 1 ' 2 ^ is the Hamiltonian of two uncoupled lay- 



ers (L = 1,2) [see Eq. (78)], the term with t± describes 



counts for adsorbates binding to carbons in random po- 
sitions {sf } in both layers. We choose C\(C2) — A^iBi) 
to guarantee that adsorbates bind only to carbons with 
coordination number z = 3. [The transport properties 
when adatoms bind to carbons in both sublattices are 
similar to those of monolayer graphene; see Sec. |III A| 
and Fig. [9] (bottom).] 

The conductivity of bilayer g rap hene follows from eval- 

r(BLG) 
'tb 



uating the Kubo formula [Eq. (79 1] with H t h —> H. 



The KPM results (summarized in Fig. |9| resemble those 
obtained previously for monolayer graphene (Fig. [8]) , but 
with important differences. (1) The formation of the im- 
purity band leads to a conductivity minimum about twice 



the value found for monolayer graphene [cr n 



2 /h 



electronic interlayer hopping, and the third term ac- 



(per layer)]. [This fact has been predicted before by co- 
herent potential appro ximat ion calculations of disorder 
in multilayer graphene! 54 1 55 1 . See Eqs. (11) and (53) in 
Refs. [521 and [55J respectively.] (2) For a high impurity 
concentration, n a( j — 5%, the conducitvity is strongly 
suppressed before actually forming the plateau; this cu- 
rious effect is rooted in the opening of a gap in bilayer 
graphene spectrum, due to the adsorbed species, uncou- 
pling the midgap region from higher energy states (see 
Fig. 3j bottom, and Fig. [4]). In this case, we can then 
speak of a "conduction gap." 

The bottom panel in Fig. [9] compares the conductiv- 
ity of monolayer and bilayer graphene for n ac j = 5%: 
away from the plateau, as carriers have energies sim- 
ilar to or higher than the interlayer coupling, we 
expect these systems to have comparable conductivities 
(per graphene layer). Our results indeed confirm the lat- 
ter point, although we found that for a very high car- 
rier density, |n c | > 20%, the conductivity of both sys- 
tems cannot be compared reliably within our KPM ap- 
proach: increasing the carrier density up to such values 
originates carrier energies close to the Von Hove singular- 
ities, and strong (spurious) numerical oscillations in the 
KPM expansion cannot be avoided. In addition, these 
oscillations behave differently in both systems (in par- 
ticular, because bilayer graphene has four such singulari- 
ties), making any comparison difficult. This is the reason 
why we have presented the conductivity for low carrier 
densities, which also coincides with the most relevant ex- 
perimental regime. 

We finish this section by noting that vacancy-induced 
disorder leads to effects similar to those reported here, a 
fact satisfactorily explained by the model of strong short- 
range scatterers presented in Sec. |III A[ For vacancies, 
though, the strong conductivity electron-hole asymmetry 
(caused by the offset resonant peaks) will not be present. 



IV. SCATTERING IN A BIASED BILAYER 
GRAPHENE 

When V 7^ 0, electrons in a graphene bilayer are de- 
scribed by Eq. (28). In this case, the energy spectrum 
develops a Mexican hat form, as represented in Fig. [10] 
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and the spectrum opens-up a gap. When the energy of 
the electrons is lower than \ V\, the Fermi surface becomes 
a ring around the Dirac point, with an inner , fc_, and an 
outer, k+, Fermi radius in momentum spaceP^^ 

Therefore, for E < |V|, we have two degenerate states 
with different momentum values. As we show below, the 
description of scattering in these two regimes, E ^ |V|, 
is necessarily different. 



The regular eigenstates of Hamiltonian in Eq. ( 28 ) in 
polar coordinates are given by 



1 



a k J m (kr) 
TbkJ m +2{kr)e 2lf> 



(82) 



(83) 



to which corresponds the eigenvalues 

E(k) = ±^V 2 (l-e k /t ± ) 2 +e 2 , 

where = Vph 2 k 2 /t± is the energy of electrons in bilayer 
graphene for V = 0, and the coefficients a k and b k read 

hl + Vil-^/t^/E} 1 / 2 , 



a, : : \i - [J - ML-- <,n_)jty ' . (84) 
h = ^-{l-V(l-e k /t ± )/E}^. (85) 
Additionally, the relation a 2 }) 2 . = e 2 k /(AE 2 ) holds. 




Figure 10: (Color online) Energy spectrum of a biased 
graphene bilayer. Several quantities denned in the text are 
depicted, and Ef stands for the Fermi energy. Information 
on the two regimes Ef ^ |V| is included. Full circles repre- 
sent degenerate states with energy E = E(k+) = E(k-), a 
fact that will have to be taken into account when establishing 
a scattering theory. 



The density probability flux Je is given by Eq. (31 1 
plus an additional term j[ , reading 



where the operator jj is given by 



fv _ 

Jo — 



-d e 
d. 



(86) 



(87) 



Throughout, we consider that electronic carriers have 
positive energy E > (the other case follows immedi- 
ately). Let us establish here some useful relations for 
later use. The energy gap A g is determined by 



A g = 2E(k min ) = 2\V\t ± (V 2 +tl) 



-1/2 



where fc m j n is defined in Eq. (|90j). Given a state with 
energy E, the two momentum values are obtained from 
the inversion of the energy spectrum, Eq. (83 1, and are 
given by the positive roots of the equation 



t± 



1 ± f(E)] 



(89) 



with f(E) = y/1 - (1 + t 2 JV 2 ){l - E 2 /V 2 ). From 
Eq. ([89]) we see that for E < \V\ the two roots are 
real, corresponding to two propagating states, whereas 
for E > \V\, only one root is real, corresponding to a 
single propagating state; this is consistent with the dis- 
persion depicted in Fig. 10 In the latter regime, the 



imaginary root is essential to fulfill the scattering bound- 
ary conditions, as in the case discussed in Sec. |II B| For 
energy E = |V|, we are at the boundary between the two 
regimes introduced above: E ^ |V|. In this case, the 
scattering descriptions below and above E — \V\ must 
provide the same answer. For E = \V\ we have fc_ = 
and k + = A 9 /(v / 2vf^) = V^mm! for E < \V\ we have 
a simple relation between fc_ and fc + , reading 



fc_ = ,2k 2 



and k„ 



A, 



2vf^- 



(90) 



The radial velocity of the electrons at fc_ and k + is given 

by 

. 2v 2 F hV 2 f(E) . , , . 

«r(*±) = -^ f^(±k ± ). (91) 

Clearly, the state with momentum fc_ has a negative ve- 
locity; the scattering formalism has to take this aspect 
into account. 

Because the regimes E > \V\ and E < \V\ are distinct, 
in the sense that the latter case contains two degener- 
ate propagating states, we develop the scattering theory 
separately for both cases. 



A. The E > \V\ regime 

For E > \V\, the two momenta are k + = k and 



fc_ = iy k^_ — 2/Cj^jjj = in. The latter value originates 
an evanescent wave at the boundary of the potential. As 
in the case in Sec. |IIB| it is simple to show that a wave 
function of the form 




1G 



represents an incoming plane wave of momentum k; = 
(fc+,0) = (k, 0) and a scattered cylindrical wave of mo- 
mentum = fc + (cos0,sin0). Note that relative to the 



case of the unbiased bilayer case, Eq. (98 1 differs in the 



presence of the amplitudes a k and b k . The scattered ra- 
dial flux has the usual form J r = v r (k)\f(9)\ 2 /r, from 
which the differen tial c ross section follows as a(0) = 
\f(6)\ 2 . As in Sec. 



HID 



we seek a wave function in the 
form of a superposition of Bessel functions of different 
kinds, which in the present case can be written as 



+ A' 



-b k e 2l6 Y m+2 {kr 



+ A' 



The ratio A 2 m j A™ reads 



,2H) 



a k J m (kr) 
-b k e 2ld J m+2 {kr) 

a k Y rn (kr 
f m +2(kr 

a K K m (Kr) 
(nr) 



K 



m+2 



(93) 



Af _ a k b lK J m {kR)K m+2 (nR) - b k a iK J m+2 (kR)K m (KR) 

AT ~ b k a iK K m ( K R)Y m+2 (kR) K m+2 ( K R)Y m (kR) 

(94) 

Combining Eqs. ( |66| and ( [94| , the equation for the phase 
shift S m follows at once. Indeed, the expression for the 
dc conductivity of electrons with Fermi momentum k + is 
similar to Eq. ( 56 ) , reading 



4e 2 



h 4n;A(fc 4 



(95) 



In the regime k + 3> v2fcmin, we have k w fc+ = k, a k ft 
ai re , an d b k ft bi K , and therefore the phase shifts given 



by Eq. (68 1 and (94 1 are essentially identical; that is, we 
have 



So 



7T 

2 



(fc+ > V2k min ) . 



(96) 



carriers (holes) must equal the phase shifts for positive 
energy carriers (electrons) if the sign of V is reversed. 

The dependence of the phase shifts on the gates volt- 
age (that is, on both k and V") is now more involved. 
Figure [TT] shows the non-zero phase shifts for electrons 
for the particular case of weak interlayer potential V. 
Similarly to the unbiased bilayer (V — 0) there are four 
(non-zero) phase shifts, however, as stressed above, the 
presence of the interlayer potential lifts the degeneracy 
observed in Fig. [6j in particular, for \V\ > the phase 
shifts with m = — 1 and m = —3 differ very much (ex- 
cept for energies very close to V). On the contrary, the 
phase shifts So and <5_2 just differ significantly close to 
the vicinity of E — V, where the systems approaches the 
"Mexican hat." 
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Figure 11: (Color online) Dependence of nonzero phase shifts 
S m on E, for the biased bilayer graphene with R — ao for low 
electrostatic potential V = 4xl0~ 3 t±. The energy range here 
excludes the interval [4, 5[ (10 -3 ix) for which the energy be- 
gins to fall within the "Mexican hat" (Fig. 10 1. In the vicinity 
of E = | V\ , we have k+ — > v2femin and the s wave phase shift 
So for electrons (or 8-2 for holes) drops quickly to the value 
indicated in Eq. ( |98| . 



As a consequence of Eq. (|96| , the conductivity is essen- 
tially linear in V g at a high electronic density. 

When the gate voltage is reduced, bringing the Fermi 
energy close to V, we have k — > 0, but k + > ^/2k m \ n is 
finite. In this case, we have 



J m+2 (kR) 



A!-, 



AT Y m+2 (kR) 



(97) 



and considering that kR < 1 , the s wave phase shift tends 
to 



g (k m i n R) 



for fcj 



v2fc m ir 



(98) 



The bias potential acts differently on electron and hole 
carriers [see Eq. (28)], with the effect that the symme- 
try relation between phase shifts changes to S m (E, V) = 
S-m-2(—E, V). Also, the phase shifts for negative energy 



B. The E < \V\ regime 

As discussed at the beginning of Sec. |IV| in the case 
E < \V\ there are two degenerate propagating states, 
characterized by k— and k + . Thus, the matrix element 
of the potential between these two states is finite, and an 
incoming particle with a well-defined momentum or 
k + ) will be scattered in a superposition of both momenta. 
This fact requires the modification of the scattering for- 
malism introduced above. 

In what follows, we develop the scattering formalism 
assuming that the incoming electron has momentum fc+; 
the case where the incoming electron has momentum fc_ 
follows immediately, and only the final results are given. 

We start by assuming that the total wave function in 



17 



the presence of the potential, at large distances from it, 
has the asymptotic form 



*(r) 



1 



at. 



Jk + x 



1 



Ofe, 



a k _ 

b k ^e™ 



■+ 
'+ 

— ik— r 



b k+ e 2 ' 1 



f ++(*)■ 



a ik + r 



(99) 



where f ++ {9) represents the scattering amplitude consid- 
ering that the outgoing electron has the same momentum, 
k + , as the incoming one, and /_| (6) represents the scat- 
tering amplitude considering that the outgoing electron 
changed its momentum to fc_. Let us stress again that 
E(k-) — E(k+). Since the velocity of the state with mo- 
mentum fc_ is negative, the sign of the argument in the 
exponential of associated cylindrical wave function has to 
be negative, since these states represent particles prop- 
agating backward in time (a positive sign gives a radial 
incoming flux). The fluxes associated with the first, sec- 
ond, and third terms on the right-hand side of Eq. (99) 
read 



J+ = v x (k + ) . 



(100) 



where 5 m ,++ is the phase shift of the partial wave m, < 
r lm.++ < 1 is a real number accounting for the transfer of 
probability flux to the outgoing momentum channel 

and < |?7 m ,H | 2 < 1. Conservation of the radial flux 

for each partial wave to imposes 



V 2 m ,+++)Vm,+-\ 2 = l- _ (105) 

Summing over to, according to Eq. (11 1, we obtain ^(r) 
in the form given by Eq. ( 99 1 , with the scattering ampli- 
tudes defined as 



/ ++ = ^=E(w + ^"--l) e '"«(106) 



(107) 



As in Sec. |II B| we write the exact partial wave of the 
full scattering problem, for r > R, as 



Jt =v r {k + )\f ++ {9)\ 2 r-\ 



and 



-v r (k-)\U-(0)\ 2 r- 



(101) 



(102) 



respectively, from which follows the existence of two scat- 
tering cross sections, defined as 

« ++ {0) = \U+(0)\ 2 and a + _(9) = - Vr ^\\f + _(e)\ 2 . 

v r (k + ) 

(103) 

Both cross sections must enter in the relaxation time 
needed to compute the dc conductivity. 

We now assume that a partial wave in the angular mo- 
mentum basis of the total wave function has, at large 
distances from the potential, the form 



a fc + 



-i{k+r— A m — m6) 



»7m,++ e 
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(104) 
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(108) 



Expanding Eq. (108) for large r and comparing it with 



Eq. (1041, we see that 



Am 

AT 



Af 
AT 



(109) 



Calculation of the differential cross section requires the 

determination of r] m ^ ++1 rj m ^ , and $ m . ++ . In the limit 

Vq — > oo, the boundary condition is , 3/ m (r = R) = 0, 
leading to 
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H^\k + R)H {2) +2 (k^R) 



b k+ a k _ H<£ +2 {k + R)H%> (k-R) 
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b k+ a k _H^(k^R)H^ 2 (k + R) - a k+ b k _H^ 2 (k_R)HX>(k + R) 



rC 1 ) 
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H^ +2 (k + R)H^'(k + R) - H^>(k + R)H^ 2 (k + R) 
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(110) 
(111) 



Although not immediately obvious, the parameters 

7 ?m,++ and J?™, h j as given by Eqs. (110 I an d ( 111 I, obey 

the flux conservation relation in Eq. (105 1. When the 



Fermi energy, Ep, approaches the energy E — \V\ from 
below, we have fc_ — > 0. In this limit, we find 



2i8 Q 



. + + 



H { 2 \k + R) 
H { 2 \k + R) 



r) m ,+- -> 0, 



(112) 
(113) 



as it should. Since fc + i? < 1, it follows from Eq. (112) 
that 



So.. 



(114) 



which, for <5o,++, gives the same result found in Eq. (98 1 



The above results hold for an incoming electron with 
momentum k + ; when the electron has momentum fc_ we 
have the cross sections: 

<r—(9) = If— (0)\ 2 and a_ + (9) 



Vr{kj, 

v r (k- 



-|/-+(0)f 



(115) 

whose amplitudes are given by the right-hand side of 
Eqs. (106 1 and (107), respectively, upon interchanging 



k + with fc_ 



C. dc conductivity of a biased bilayer graphene 



As discussed in Scc. |IIID| calculation of the dc conduc- 
tivity requires the computation of the exact phase shifts. 
We start by studying the behavior of the s wave phase 
shift as a function of the Fermi momentum for a biased 
graphene bilayer. 

In the biased bilayer, the ability of independently tun- 
ing the electronic density and the value of the gap A g 
requires the use of two gates, a bottom and a top gates, 
as shown in Fig. [12] The electric field in the top-gate 
dielectric is (e > 0) 



E, = 



en t 
eteo 



(116) 



and that in the bottom-gate dielectric is 



£b = — , (H7) 
where n t and rib are the electronic density in the top 




Figure 12: (Color online) Capacitor geometry for dual-gate 
transistor.^ The figure is self-explanatory. Values of the sev- 
eral quantities are: 8 = 3.4 A, b — 300 nm, and t = 20 nm. 
V t and Vb stand for the top and bottom gate potentials, re- 
spectively. 



and bottom gate, respectively, and et and eb are the rel- 
ative permittivity of the top- and bottom-gate dielec- 
tric, respectively. Charge neutrality requires that the 
total amount of charge accumulated in the bilayer is 
—en = — e(nt + Ht>)- The electrostatic potential differ- 
ence between the top gate and the bilayer is V t = tE t , 
whereas between the bottom gate and the bilayer it is 
V b = bE h . It follows from Eqs. (flTol) and (111 71 that 



V h =b 



£b£o 



ben 
£b£o 



bet 



V t 



(118) 



bilayer is given by 



Inverting Eq. (118 1, the total electronic density in the 



n = V h 



be 



et 



Vt. 



(H9) 



When n is positive, the bilayer is doped with electrons; 
when n is negative the system is doped with holes. Fi- 
nally, the electrostatic potential difference between the 
two graphene layers in the bilayer is given by 



AV = (E h - E t )S 



neS 
£b£o 



ft 



+ 1 )-Vt , (120) 



where 6 = 3.4 A is t he inter layer distance (we are ignor- 
ing screening effects, ^ * 37 * 39 ! which are not important for 
small V t ). The variable V introduced in Eq. (26) relates 
to AV as 2V = AV. Taking typical values for dual-gate 
bilayer transistorspS we have: e Si o 2 = 3.9, £Hf0 2 = 25, 
£nfc = 2.4, b = 300 nm, and t = 20 nm (both dielectrics, 
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Hf02 and NFC, have about the same width). The rela- 
tive permittivity of e t is 



ft 



2eHf0 2 e NFC 



£Hf0 2 



£nfc 



(121) 



In working devices p^J we have |Vb| < 70 V and |Vt| < 4 
V. 

The calculation of dc conductivity follows, as before, 
from Boltzmann's transport theory. In the regime E > 



|V|, (Tdc is still given by Eq. (56 1, but with the phase 
shifts determined from Eq. (94 1. When E < \V\, there 
are two scattering channels and this implies that the re- 
sulting formula for o^c differs somewhat from that given 



in Eq. (56 I, reading 
4e 2 1 



± 1 

rii<j(k + ) n,cr(fc_) 



(122) 



h 2 

where a(k±) is defined as 



/•27T 

a{k±) = [cr ±)+ (0) +<t ± _(0)](1- cos 0). (123) 
Jo 

Inserting the expressions for the differential cross sections 
[Eqs. (1031 and (115)] and performing the integral yields 



r(fc± )^Re£[j 



Vm,±±e 



2iS„ 



1 - (?7m,±±e 



2i6 n 



1) (Vn 



■+l,±± e 



-2iS m +i,±± 



1) 



k± ty(fczp) , 2 _ * \ 



(124) 



The formulas for cr and a |_ are identical and thus 

are not presented. The dc conductivity follows from the 
determination of the Fermi momentum, given the carrier 
density in the bilayer ; wh ich in turn depends on both 
gates as given by Eq. (119 1. The relation kp = irn (valid 



for various two-dimensional systems) must be adapted 
to take into account the degeneracy of the spectrum 
(Fig. 10 1 and reads, 



-n + k min , 



(125) 



and the other propagating state k F relates to kp accord- 
ing to Eq_j90|. 

Figure |13| shows the dc conductivity as function of the 
back gate for fixed values of Vt. As the back gate Vb is 
varied, the gap A g and the Fermi energy change; for a 
small window of width ~ IV around Vb — —17 V the 
system moves into the regime E < |V| and expression 
in Eq. ( 125 1 must be used to determine the carriers en- 
ergy. In this energy regime, kp is bounded according to 
\/2fcmin > kp > fe m m and hence the value kp — is for- 
bidden; as a consequence, and at odds with the unbiased 
bilayer, the minimum conductivity is not exactly zero, 
having a value of <J m i n ~ 3e 2 jh for Vt = 1 V. 



V. CONCLUSION 

In the early studies of transport in graphene, charged 
impurities located in the substrate seemed to explain 
the measured conductivity. Recent experiments suggest 
other possibility thoughJ2212il While there is a consen- 
sus that electron and hole puddles, induced by charged 
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Figure 13: (Color online) Dependence of the biased bilayer dc 
conductivity on the back-gate potential Vb (values of Vt are 
indicated). The solid line shows the dc conductivity for the 
unbiased case V = 0, for comparison. 



impurities, dominate the landscape near the neutrality 
point, away from this point, adsorbed hydrocarbons, at 
the surface of graphene, may be the limiting factor in 
dc-transport. 

In the present paper, we established an intuitive theo- 
retical picture of scattering due to resonant scattering 
originated by adatoms. Although resonant scatterers 
have been studied before (first in Refs. 1411421 and, more 
recently, in Ref. l20)l . we have established the first coher- 
ent picture of resonant-scattering limited dc transport 
valid for both monolayer and bilayer graphene. 
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Section IIII Al reviews the electronic structure of mono- 
layer graphene and presents, for the first time, the den- 
sity of states of bilayer graphene with resonant contam- 
inants. Despite the distinct electronic structure of pris- 
tine monolayer and bilayer graphene, this section shows 
that resonant adatoms lead to the same effect in both 
systems: the emergence of resonant peaks in the vicinity 
of the Dirac poi nt, a situation reminiscent of vacancy- 
induced disorder ! 26 l 44 l Using a simple tight-binding toy 
model, resonant adatoms are seen to be reliably mim- 
icked by a particular class of short-range scatterers, that 
is, those having an intrinsic energy much higher than 
typical graphene energies. This fact motivates the sub- 
sequent study of dc transport using strong s hort-ra nge 
potentials in a continuum formulation (Sees. 

TTTd i. 



IIIC and 



Section |III| shows that the typical dependence of con- 
ductivity with the electronic density in the monolayer 
(sublinear dependence) and bilayer (linear dependence) 
systems can be explained assuming resonant scatterers 
alone. The comparison with experimental data bears out 
the agreement with dc transport experiments performed 
in exfoliated few-layer graphene films, hence providing 
further strength to the resonant-scatterer hypothesis. To 
justify the robustness of a continuum-model approach 
based on strong short-range scatterers as prototypes of 
real resonant adsorbates, we have calculated the semi- 
classical conductivity due to two types of strong local 
potentials (hard-disk and delta-potential) , finding perfect 
agreement between the two methods (partial wave analy- 
sis and Lippmann-Schwinger equation, respectively) and 
tested the validity of the long wavelength limit (on the 
basis of the continuum formulation) against numerical 
lattice calculation using a T matrix approach (Sec. Ill El. 



Section IIII Fl demonstrates the incorrectness of the 
widely used FBA within the semiclassical (Boltzmann) 
approach, in the context of short-range disorder, and the 
need to compute the electronic scattering amplitudes as 
accurately as possible, hence, clarifying an issue over- 
looked in the graphene literature. Section fill G| presents 
the Kubo dc conductivity evaluated numerically with a 
KPM; from this calcution, the breakdown of the semiclas- 
sical picture close to neutrality, in the regime of a high 
concentraton of impurities, is clearly observed. Here, the 
case of bilayer graphene is addressed for the first time, 
with the results showing that a "conduction gap" takes 
place for selective adsorbate bonding, due to a strong 
supression of the conducitivy in the surroundings of the 
resonant impurity band. 

Finally, due to its importance for technological appli- 
cations, scattering in the bilayer graphene with a gap 
in the sectrum is studied in Sec. |IV| by extending the 
well-established partial wave method (Sec. [TTJ) to describe 
scattering in the biased bilayer graphene. Such a scatter- 
ing theory has never been developed before (to the best of 
our knowledge) and can be easily adapted to tackle other 
physical scenarios requiring the need for computing scat- 
tering amplitudes when the energy dispersion relation is 



degenerate. 

We are confident that our results help to elucidate the 
electronic transport properties of this remarkable two- 
dimensional material. 

Note added: After submission of this work for publica- 
tion, we become aware of a paper 77 which also discusses 
the effect of resonant scatterers on the dc conductivity of 
single-layer and bilayer graphene, with results that are 
consistent with ours. 
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VI. APPENDIX 

In this appendix, starting from the low-energy contin- 
uum theory, we derive the nonperturbative semiclassical 
dc conductivity of monolayer and bilayer graphene with 
short-range scatterers. This calculation requires the so- 
lution of the two-dimensional scattering problem, where 
a massless fermion with incident momentum p = Kk is 
brought to interaction with an impurity. We model 
the potential of the impurity by a delta function, V c i — 
Vq8(y). Following standard methods, the formal solution 
of (Hq + Vd — -E^k = can be written as, 



*k = 0k + GoV d V k . 



(126) 



where Hq is the low-energy Hamiltonian of graphene; 0k 
is the solution of the free problem (Ho — E)(pk = and 
describes the state of the incident particles. Here, Hq 
refers to the Hamiltonian obtained from expansion of the 
graphene dispersion around the K point (the calculation 
involving the remaining valley is equivalent). The resol- 
vent is given by Go = 1/(E + i0 + — Hq), and the energy 
includes a small positive imaginary part i0 + . The spinor 
</>k(r) = (r|0k) has the forrrP, 



0k (r) 



with. 



'2A 



1 



(127) 



(128) 



The Berry phase is ifB = 7rA and equals 7r for mono- 
layer graphene, whereas for bilayer graphene its value is 
2-7T [compare with Eq. The second component of u k 
includes the sign s — ± of the electronic carrier charge 
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and 0k = arctan(fc 1/ /fc a; ). Switching Eq. ( 126 I to the posi- 
tion representation, we obtain the Lippmann-Schwingcr 
equation, 

*k(r) - 0k(r) + j d 2 v'G {v - r')V(r')*k(r') • (129) 

In the latter equation, Go(r — r') = 
(r| (E + iO+ -Hof 1 |r') is the Green function of 
the problem. Monolayer graphene has Hq — hvper ■ p 
and the Fourier transform of the Green function obeys, 



(E + i0 + - ct • p) Go(p) = 1. 



(130) 



where Go(p) = / dv exp (— ip ■ r) Go(r) (notice that to 
simplify notation, we have set h = 1 and Vp — 1). In- 
verting the 2x2 matrix on the left-hand side of Eq. ( 130 I , 
we arrive at 



G (p) = ff i(p)(£ + <x-p) 
ffi(p) = l/[E 2 -p 2 + iO+ 



(131) 
(132) 



The calculations for E > and E 1 < are similar and 
to be specific we focus on the former situation. Indeed, 
the inclusion of a small imaginary part from positive val- 
ues i0 + amounts to consider outgoing waves (see below). 
We write E = k and evaluate the Green function in the 
real space representation, 

G (r - r') = JL (£ - Ur ■ V) f dV^* -0 5l (p) 

(133) 

= -i(fc-« T .V)ff (1) (k|r-r'r), (134) 

where Hn^ (k\r — r'|) is the Hankel function of the first 
kind of order n, whose asymptotic form is that of outgo- 
ing cylindrical waves [see Eq. ( 149 1]. The Hankel function 



obeys 8 x Hq(x) + H{ 1 '(x) = 0; hence, 



a • VH^(k\r - r'|) = -kH[ L \k\v - r'\)a e , (135) 
where we have introduced the matrix, 



equal to the unit; the Fourier transform of bilayer Green 
function reads 



G (p)=ff 2 (p) S + 7 cr.D(p) 



(138) 



where 7 = l/t±, E = 7 fc 2 , D(p) = [p 2 x - p 2 y , 2p x p y ) T , 
and 



92 (p) 



1 

2E 



1 



1 



E --fp 2 + i0+ E + jp 2 + i0+ 



(139) 



Since Hq is quadratic in momentum operators, 32 resem- 
bles a non relativistic propagator. Again, we focus on 
the case of electrons (7 > 0) , 



G (r -v') = ^L(k 2 -<r-T>) j dV^-'^p) . 

(140) 

The contribution to the integrand of Eq. ( 140 1 with poles 
in the real axis can be simplified using, 



I ?7T 1 

= ^[S( P + k) + S(p-k)}+P.V 



k 2 -p 2 + i0+ 2k 



Performing the integral in Eq. ( 140 ) yields 



k 2 — p 2 
(141) 



G„(r-r'i = -l^^- a .-D)\-iH^(k\r-r'\) 



-K (k\v-r'\ 

7T 



(142) 



The first term in brackets describes scattered waves in 
two dimensions, whereas the modified Bessel function K 
describes evanescent waves (recall that k — > ik is a solu- 
tion of Hq with the same energy) . For short-range poten- 
tials the main contribution to the scattering amplitude 
comes from evaluating Eq. ( 129 I within the region where 



I r — r' I 3> 1 and hence Kq will not contribute (see later) . 

In what follows, we compute the nonperturbative scat- 
tering amplitude for monolayer and bilayer graphene, 
which will be needed for the calculation of the dc con- 
ductivity in these systems. 



<T0 




(136) 



and the angle 9 = 6(r,r') is defi ned b y (r — r') /|r — r'\ = 
(cos#,sin#) T . Combining Eqs. ( 135 1 and (134), we have, 
at once, 



G (r- 



ik r 



(137) 

The derivation of the Green function of bilayer 
graphene follows identical steps. We write the free Hamil- 
tonian as H = —(vph 2 /tj_)(T ■ D, with D = (d 2 — 
d 2 , 2d x d y ) T . As before, we set h and vp temporarily 



Al. Nonperturbative amplitude for monolayer 
graphene 

Inserting the expression of the potential Vd = VoS(r) 
in the Lippmann-Schwinger equation [Eq. ( |129 1] and per- 
forming the spatial integration results in, 

*k(r) = 0k(r) + ^ o G o (r)* k (O), (143) 

which is ill defined because putting r = yields a diver- 
gence, namely, 5'k(0) = 0k (0) + (00). This stems from 
the singularity of Gq(t) [Eq. (137)] at the origin r = 0, 
a common situation in field theories. The only way of 
curing this divergence is by means of renormalization.^1 
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Let us explicitly describe this procedure. The explicit [Eq. (148)] allows us to find the explicit expression of 



expression for Gq(0) [see Eq. (133)] reads 
d 2 p 



G o (0) 



(27T)' 



gi (\p\)[E + a-p] , 



(144) 



Evaluating Go(r) at the origin and setting E — k yields 



Go(0) 



dp-. 



P 



(145) 



which is logarithmic divergent. To obtain a physical 
meaningful result, a momentum cutoff, J3 max , in the up- 
per limit of the integral must be considered. (Such pro- 
cedure is justified because graphene, being a solid, has 
an intrinsic energy cutoff of the order of the bandwidth.) 
We thus have 



G o (0) = 



k f p 
2tt 



dp 



P 

k 2 -p 2 + i0+ 



This integral yields 



G (0) = — ln(kR) 



(146) 



(147) 



where we have assumed k <C p maK and R = 1 /p max is a 



length scale of the order of ao- Setting r = in Eq. ( 143 1 
using the latter result and solving for x Pk(0) gives 



*k(0) 



Vo 
2tt 



kln(kR) 



c(0) 



1 - — klntkR) 
2tt v ' 



,(!) 



(148) 



To identify the scattered amplitude, we need the 
asymptotic form of the Lippmann-Schwinger equation 
[Eq. ( 129 )]. For short-range potentials the main contribu- 
tion in Eq. ( 129 ) comes from the region where |r— r'| 3> 1. 
Inserting the exact form of the propagator in space rep- 



resentation [Eq. (137)] and using 
H^\k\v-v'\)^ 



ik\r— r | 



FfVfclr-r'l) -> -ij : -e 

1 K 1 u Y iirk\r-r 

leads to 
* k (r) = k (r) 



ik\r— r' | 



(149) 
(150) 



ik 
8nr 



riVe"*—' (l + <7 )V d (r')* k (r'), 

(151) 



where we have approximated |r — r'| ~ r — r • r'/r and 
identified the wave vector at the point of observation, 
k out = kr/r. The exact form of the spinor at the origin 



# k (r); letting 5 e = a e {v' = 0), 
*k(r) = M r ) 

_ Vo rw t 

1 - ^k\n(kR) V 87rr' 



lkr (1 + a e ) 4 1J • (152) 



The action of (1 + a$) on the spinor u^' yields the 
Berry phase term for scattering in graphene; without loss 
of generality, we take the incident momentum along the 
x axis, k = (k, 0), and thus 



1 



X 1 ) 



where 



E B (6) = {l + e- ie ) 



(153) 



(154) 



and the scattering angle reads 9 = Z (k, k out ) [recall 
Eq. (136 1 and comments therein]. 

The wave function of the scattered particles is then 



fkr 



(1) 



f k (r)^(r)+/(9) T « ko 
with the scattering amplitude reading, 



(155) 



/(*) 



1 

hvp 



Vo 



8tt 1 



Vo 
2-jtIi v j 



-k\n(kR) 



Eb(0), (156) 



where we have restored all the constants. (Note that here 
Vo has units of [energy] x [length] 2 ; the relation between 
Vq and the effective impurity potential V e g in a lattice 
theory can be shown to be Vo ~ A c V e s, where A c is 
the area of graphene 's unit cell.) This result is to be 
compared with the result from the FBA, which amounts 
to approximate 5 , k (r') by the unperturbed wave function 
k (r') in Eq. pi): 



/Born(#) 



1 



-v E B (e). 



(157) 



The latter is only accurate in the limit of a very small Vo, 
which is of limited interest. The nonperturbative result 
discloses a singular momentum, fc s ; ng , 



^sing In (^sing^) 



2nhvF 
Vo 



(158) 



which corresponds to a bound state of our problem. More 
importantly, the nonperturbative amplitude for Vo — > oo 
[recall that resonant scatterers in graphene give origin to 
strong short-range potentials, see Sec. ( III A l] reads 



fv^ooiO) = 



E B {6) 



2 y/k\n(kR) 
which is the main result of the present section. 



(159) 
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A2. Nonperturbative amplitude for bilayer graphene 

Calculation of the scattering amplitude for the bilayer 
graphene follows as in Sec. Al, albeit with the important 
difference that the Green function does not diverge at the 
origin and thus no renormalization procedure is needed 
this time. This explains why no regularization length 
appears in the final result for the conductivity of bilayer 
graphene. We now outline the derivation of this result. 



(2) 

The action of the last term on the spinor u k yields the 
bilayer Berry phase term; taking the incident momentum 
along the x axis, k = (k, 0), we obtain 



(l + *2*)u£° = S s (2#) ' 



V2A 



1 

JliB 



ee E B (20)v™ ut , 



The explicit expression for Gq(0) [see Eq. (140)] reads 



Go(0) 



d 2 p 



ga(p) £ + 7 <T-D(p) 



(2nY 

which, setting E = jk 2 , can be simplified to 
d 2 p 1 



G o (0) - 7 fr 



(2tt) 2 ( 7 fc 2 -M0+) 2 - 7 V 



(160) 



(161) 



(167) 



where Eb is defined in Eq. ( |154[ ) and 9 is the scattering 
angle, 9 = Z(k, k out ) [recall Eq. (1361 and comments 
therein]. The wave function of the scattered particles is 
then 



*k(r) = </> u (r) + /(0W.* k 



e* kr (2 ) 



with the scattering amplitude reading 

r 2v 



The above integral can be solved straightforwardly by 
contour integration; the result is 



/(*) 



2kir 8v 2 F h 2 /t 



iVr 



S B (20) 



(168) 



(169) 



G o (0) 



i 

'87 



(162) 



where we have restored all the constants. The FBA is 
recovered in the limit Vq <C energy scales, 



The amplitude of the wave function at the origin 
[Eq. ( 143 1] therefore reads 



/Born(#) 



v 



±v 2 F h 2 /tj 



;Hb(20). 



(170) 



*k(0) 



1 + i 



.Vo 

87 



c(0). 



(163) 



l± v 2kn 

In contrast, in the limit of interest Vo — > 00, we obtain 

T~ B (2#) 



jV->oo(#) — - 



To identify the scattered amplitude, we have to repeat 
the derivation of the asymptotic form of the Lippmann- 
Schwinger equation [see Eqs. ( 149 I ( 1 65 l] . The asymp- 



«7T y/k 



(171) 



totic form of the propagator can be calculated from 
Eq. p2l, 



A3. The dc-conductivity of monolayer and bilayer 
graphene 

The dc conductivity follows from the Boltzmann equa- 



Go(r-r') 



1 i f 1 e -2i0(r,r') 

2kn\r- r'| I e 2t9 ^' r "> 1 



(164) 

where 9 — Z(r,r'). Inserting the latter expression into 
Eq. (1291, and approximating |r— r'| ~ r— r-r'/r, permits 



us to identify the wave vector at the point of observation 
3 k ou t = kv/r, 



tipn (see Sec. Ill B I . The expression for the semiclassical 
current j can be manipulated to yield a more convenient 
form of the conductivity for our purposes. We reproduce 
the main steps; at T = the Fermi function becomes the 
Heaviside function, 9{tk F ™ £fc), and hence the expression 
for the current (including spin and valley degeneracies) 
reads 



. _ g s g v e 

3 — l^-\2 



* k (r) = k (r) 



Akr 



(2n) 2 



J dk T (k)6(e kF - e k )(w k ■ E)v fc . (172) 



4 7 V 2knr 
x (l + <7 2fl(r , r0 ) Vd(r')* k (r'), 



d 2 r / e -ik out -r' 



(165) 



Performing the angular integration, and using the rela- 



where the definition of <jg is given in Eq. ( 136 1. As before 



letting (T20 = crg( r r /_ ), and using Eq. (163 1, we get 



e 2 k F 
ttH \v r (kp] 



Mk F )(v kF -E)v fcF . (173) 



* k (r) = k (r) 

87 + Wo V 2knr 



The longitudinal dc-conductivity follows from the latter 
expression: 



e lkr (l + a 2g )u^ . (166) 



2e 2 

(Tdc = -rT(k F )\v r (k F )\k F . 
h 



(174) 
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Using the results in Sees. Al and A2 and the definition of 



relaxation time r(kp) (Sec. Ill B I , we can readily obtain 



the dc conductivity in the regime of Vo ^energy scales. 
(For a discussion of the on-site energy Vo magnitude see 
Sec. [inA| ) 

The dc conductivity in the limit Vo — > oo reads 



strong 
7 dc 




In 2 (kpR) for monolayer 
for bilayer 



(175) 



As expected, the dependence on kp coincides with that 
obtained through the partial wave expansion method em- 
ployed in Sec. |III| Remarkably, the expressions match ex- 
actly [compare with Eq. (64) and (72)]. This entails that 
scattering of a hard disk of radius R ~ ao and scattering 
off a strong delta potential have the same dependence on 
the momentum of the incident particles (in both mono- 
layer and bilayer graphene). 
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